##################################
########Models t-1 to t-10########
##################################

#Figure A.1 

#Full Model ALL
library(sandwich)
library(interplot)
library(AER)
library(DAMisc)
library(multiwayvcov)
library(tidyverse)
library(ggplot2)
library(wrapr)
library(jtools)

##############
###FIGURE 1###
##############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) +
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1","lag.lngcp",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction")])

omit<-NDs
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("Figure1.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist",
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


################
###FIGURE A.1###
################

#############	
#############	
### t – 1 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1","lag.lngcp",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc1.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 2 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 2)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction2")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction2+lag.lncapdist+lag.cost_sanction2:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc2.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction2", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 3 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects(t - 3)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lnNTL","lag.lnbdist1","lag.lngcp","lag.lnpop_gpw_sum",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction3")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction3+lag.lncapdist+lag.cost_sanction3:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc3.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction3", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 4 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 4)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction4")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction4+lag.lncapdist+lag.cost_sanction4:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc4.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction4", moderator = "lag.lncapdist",
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 5 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 5)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction5")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction5+lag.lncapdist+lag.cost_sanction5:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc5.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction5", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 6 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 6)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction6")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction6+lag.lncapdist+lag.cost_sanction6:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc6.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction6", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 7 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 7)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction7")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction7+lag.lncapdist+lag.cost_sanction7:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc7.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction7", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 8 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 8)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction8")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction8+lag.lncapdist+lag.cost_sanction8:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc8.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction8", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 9 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 9)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction9")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction9+lag.lncapdist+lag.cost_sanction9:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc9.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction9", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 10 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 10)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL", "lag.lngcp","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.cost_sanction10")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction10+lag.lncapdist+lag.cost_sanction10:lag.lncapdist+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lngcp+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("fullsanc10.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction10", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

